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Abstract 

Reactive molecular dynamics simulations are used to study initial stage of pyrolysis of ablation 
materials and their composites with carbon nanotubes and carbon fibers. The products formed during 
pyrolysis are characterized and water is found as the primary product in all cases. The water formation 
mechanisms are analyzed and the value of the activation energy for water formation is estimated. A 
detailed study on graphitic precursor formation reveals the presence of two temperature zones. In the 
lower temperature zone (<2000 K) polymerization occurs resulting in formation of large, stable 
graphitic precursors, and in the high temperature zone (>2000 K) polymer scission results in formation 
of short polymer chains/molecules. Simulations performed in the high temperature zone on the phenolic 
resin composites (with carbon nanotubes and carbon fibers) shows that the presence of interfaces had no 
substantial effect on the chain scission rate or the activation energy value for water formation. 



1 . Introduction 


Ablative materials are thermal insulators used in hypersonic space vehicles which absorb heat in part 
through endothermic, pyrolysis reactions. The char produced as a result of these reactive processes 
produce a thermally insulating and protective layer at the ablative material surface. Mechanical strength 
of the char layer is required to prevent premature flaking and subsequent exposure of the virgin ablative 
material underneath. The amount of ablative material should be enough to compensate for any part of 
the char layer that might be ablated away during the re-entry. However, the weight of the ablative 
material should be kept to a minimum. Optimization of the thermal protection system requires accurate 
prediction of the pyrolysis process and the evolution of char morphology. 

One recently developed ablative material is “phenolic impregnated carbon ablator” (PICA) which is a 
carbon liber composite with a phenolic resin matrix. This low-density material is a good choice due to 
its low thermal conductivity and efficient ablation properties. Phenolic resins have a high char yield 
under pyrolytic conditions. These reactions involve a complex sequence of reactions, and further 
development of these materials will require fundamental understanding of these molecular processes, 
and subsequent char evolution under transient conditions. Computational modeling of the resin-to-char 
process using atomistic-level simulation technique such as Molecular Dynamics (MD) with reactive 
force-field can provide insights into the involved reaction barriers and pathways. 

Recently, De-en Jiang et.al., [1] performed MD simulations with the reactive force field, ReaxFF, to 
simulate the initial stages of phenolic resin pyrolysis. Results from these simulations are benchmarks for 
the studies presented here. The time scale accessible for routine MD simulations is limited to the 
nanosecond range, which allows study of only the initial stage of pyrolysis. Lawson and Srivastava [2] 
developed a high temperature, pyrolytic MD simulation method for polyethylene, based on the gradual 
removal of free Hydrogen atoms from the system. With their method, they were able to run the reactions 
to completion, however because of the approximations involved, accurate information about reaction 


kinetics were difficult to obtain. 



In an effort to improve the mechanical properties of the char layer, NASA Ames Research Center is 
recently explored the effect of adding of carbon nano tubes (CNT) to PICA. Tensile test results showed a 
35% increase in strength with the addition of 0.4 wt% functionalized CNTs in PICA, without any 
change in insulation properties [3]. This increase in strength at low concentrations is possible because 
the percolation threshold for CNT composites is very low, on the order of 0.1 vol %. The CNT-polymer 
interface is expected to play a major role in maintaining the insulation properties by limiting the 
effective conductivity of CNTs in the composite due to a large interfacial thermal resistance [4], Thus, 
without an appreciable increase in the weight or thermal conductivity, CNT addition can provide large 
mechanical reinforcement to the composite. In addition, the char morphology is known to change due to 
the addition of CNTs, as the char strongly adheres to the CNT yielding cigar- like structure [5]. 

Furthermore, addition of CNT in flame retardant materials has been recently explored. Research has 
shown that the addition of only 0.5 wt% concentration of single-walled CNT (SWCNT) resulted in 
formation of a crack-free network-structured protective residual layer in PMMA composites and also 
reduced the heat release rate [6]. Formation of this crack- free protective layer is essential to prevent 
exposure of virgin ablative material underneath. However, in both the cases (ablation and flame 
retardant), the role of CNT in altering the (already) complex pyrolysis reaction is not known. 

In this paper, we attempt to understand the role of CNT and carbon fibers on the initial stage of 
pyrolysis. We performed detailed simulations with the ReaxFF method to characterize the products 
involved and the kinetics of pyrolysis. For the pristine phenolic resin case, (i.e., without carbon fillers) 
we also performed simulations to determine the conditions and processes involved in the formation of 
graphitic precursors. Two temperature zones were found to be important: in the lower temperature zone 
(<2000 K) polymerization occurred resulting in formation of stable graphitic precursors, and in the high 
temperature zone (>2000 K) disintegration resulted in formation of short polymer chains/molecules. The 
effect of two fillers, carbon fibers and CNT, on the initial stage of pyrolysis was also studied. No 
substantial effect was found due to filler addition on the product (water) formation rate or polymer 
disintegration rate in the high temperature zone. The paper is organized as follows: Section 2 describes 



the simulation method and details; results are presented in section 3; and section 4 provides a summary 
and conclusion of our work. 

2. Method 

MD simulations were performed with the LAMMPS molecular dynamics code [7] . The reactive force 
field, ReaxFF [8], was employed to perform pyrolysis simulations. ReaxFF uses an empirical 
relationship between bond orders and bond lengths to allow smooth bond formation and breaking. This 
method has been quite successfully applied to a variety of chemical problems, including thermal 
degradation of polymers [9]. Simulations were performed on a non cross-linked phenolic formaldehyde 
resin (figure 1: top panel) in which two, neighboring phenol moeities are connected by a methylene 
group. The simulation cell, shown in figure 1: bottom panel, consisted of 16 phenol formaldehyde 
chains with 8 repeating units each in the ortho-ortho sequence and a methyl group termination. The 16 
polymer chains were randomly placed in the simulation cell by using the software Xenoview [10], 
which uses a Monte Carlo algorithm to build the amorphous structure. The simulation box size was 
adjusted to 2.62 nm x 2.62 nm x 2.62 nm to yield the experimental resin density of 1 .25 gm/cc at 300 K. 
The polymer, simulation cell size and density are similar to the work performed by De-en Jiang et.al [1]. 
The time step used to perform the simulation was 0.25 femtosecond (fs). Periodic boundary conditions 
were applied in all the directions of the simulation cell. A Berendsen thermostat [11] with a temperature 
damping constant of 100 fs was used to maintain constant temperature of the simulation cell. The 
ReaxFF parameters for carbon, hydrogen, and oxygen were taken from a previously optimized data set 
by K. Chenoweth et.al. [12]. 

For composite systems, we chose to model two fillers - CNT and carbon fibers. Two CNT composite 
systems with two different single-walled CNT chirality, (5,5) and (10,10) were studied. The radius of 
(5,5) and (10,10) CNT is 3.5 A and 7.0 A, respectively. In both systems, the CNT is 27 A long and 
uncapped (or open) as shown in the inset in figure 2a. The phenolic resin structure is similar to the 
pristine case. The phenolic resin chains were introduced randomly around the CNT. Periodic boundary 



conditions were applied in all the directions of the simulation cell. The size of the simulation cell was 
adjusted such that the resin density was 1.25 gm/cc. The excluded volume due to the van der Waal’s 
distance between the CNT and resin was also accounted during this adjustment. Our simulations 
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indicate that the van der Waal’s distance between the resin and CNT is 3.4 A, similar to the work 
published by Wei and Srivastava for CNT -poly ethylene composite modeled with a Tersoff-Brenner 
potential and united atom model potential [13] . 

Currently, the ablation materials employed on space vehicle are made of phenolic resin impregnated 
within carbon fibers. Owing to the nanometer length- scale of the present study, the carbon fibers are 
simulated as surfaces with zero curvature, i.e., planar graphite sheets. As shown in the snapshot in figure 
2b, the graphite block has 3 graphene sheets in the A-B-A stacking sequence. The surface area parallel 
to the sheet has the dimensions, 25.56x24.6 A . The value of the surface area for the graphite block was 
selected to match the surface area of CNT (10,10). This was intended to provide insights on the effects 
introduced by the presence of tube curvature. The phenolic resin structure was similar to the pristine and 
the CNT case. The phenolic resin chains were introduced randomly on both sides of the graphite sheet. 
Periodic boundary conditions were applied in all the directions of the simulation cell and the size of the 
simulation cell was adjusted such that the resin density was 1.25 gm/cc. 

The equilibration procedure for all systems (pristine and composite) involve (i) equilibration at 0.1 K 
with an NVT-MD simulation for 40 picoseconds (ps), followed by, (ii) thermal annealing procedure - 
heating the system with a ramp rate of 0.02 K/fs to 1000 K and then equilibrating at 1000 K for 20 ps 
and then cooling to 300 K with the same ramp rate. This annealing procedure was performed twice 
followed by an equilibration at 2000 K for 200 ps. Our in-house code was used to examine and identify 
the new molecules that were produced during the pyrolysis simulations. Our analysis confirmed that no 
reactions took place during the 200 ps equilibration run. Pyrolysis simulations were performed at 5 
different temperatures - 2750, 2875, 3000, 3125 and 3250 K to obtain reaction rates and activation 
energy. To improve the statistical accuracy and provide error bars for the reaction rate values, all 
simulations were performed with 4 different starting structures at each temperature. 



3. Results and Discussions 


This section is divided into two sub-sections - pristine phenolic resin and phenolic resin composites. 

3.1 Pristine phenolic resin 

Under this subsection, first the product evolution is discussed including the associated reaction rates 
and activation energies. Then, the details on graphitic precursor formation are presented. 

A. Evolution of Products 

A detailed structural analysis on the equilibrated phenol formaldehyde resin is performed to analyze 
the accuracy of the structure yielded by the relatively new ReaxFF potential. NPT simulation performed 
to estimate the resin density at 300 K indicates that the resin density is only 4.5% lower than the 
experimental resin density of the commercially available Bakelite novolac. However, in all the 
simulations the polymeric resin density was maintained at 1 .25 gm/cc to directly compare the pyrolysis 
results with Jiang’s [1] work. Schurmann and Vogel [14] performed simulations on a similar phenol 
formaldehyde resin containing 8 repeat units with 57% in ortho-para and 43% in ortho-ortho sequence 
at 300 K and resin density of 1 .25 gm/cc with the forcefield - pcff . The radial distribution function (rdf) 
plots at 300 K averaged over a 40 ps trajectory from the amorphous phenol formaldehyde resin structure 
generated by ReaxFF (figure 3) matched qualitatively with the rdf plot in figure 3 of Ref. 14. They both 
show the length of the short-range order (due to the amorphous nature of the structure) as 0.6 nm. 
Quantitatively, compared to the polymer consistent force field (pcff), the peaks indicating the C-C, C-H, 
C-0 and O-H bonds are higher and less broad with the ReaxFF potential,. Further estimation of the 
fidelity of the ReaxFF method to simulate phenol formaldehyde resin is performed by comparing the rdf 
from ReaxFF with ab-initio MD simulations at 300 K. They show good agreement validating the 
ReaxFF potential. The details on this comparison are provided in Appendix A. 

High temperature simulations at 5 different temperatures in the range of 2750 to 3250 K are 
performed on 4 different equilibrated structures. Water is the primary product in all these simulations. 



The second most abundant product is hydrogen gas. The number of water molecules formed (averaged 
over 4 different starting structures) is plotted as a function of time in figure 4a. The rate of water 
formation (i.e., the slope) increases with temperature, suggesting a thermally activated reaction process. 
At each temperature, the value of the water formation rate is calculated from a least-square fit of the 
data sets and plotted in figure 4b. The error bars in fig. 4b indicate the difference between the maximum 
(minimum) rate from the different starting structures and the average rate. Using the Arrhenius equation, 
a value of 124 ± 20 kJ/mol is obtained for the activation energy for water formation. This value matches 
well with the value of 144 ± 28 kJ/mol published by Jiang et.al. [1], 

Additionally, simulations are performed for an equilibrated structure at three different temperatures, 
2750, 3000 and 3250 K with simulation time extended to 200 ps to reach completion of reactions. As 
shown in figure 5, initially the water formation rate depends on the temperature. At all temperatures, the 
total number of water molecules formed reaches a plateau after a certain time and converges to a value 
of ~100. The simulation cell contains 128 oxygen atoms, which indicates a yield (defined as the ratio of 
actual product formed to the theoretical) of ~78%. The different mechanisms involved in water 
formation (including the transition states) are analyzed and presented in Appendix B. The remaining 
oxygen atoms are consumed during formation of CO, C0 2 , and, aldehydes, ketones and alcohols 
(produced in small amounts). 

B. Graphitic Precursor Formation 

The ablation process during the re-entry of the space vehicles involves carbonization of the phenolic 
resin under an inert atmosphere at high temperatures ranging from 1500 to 2000 K that leads to 
formation of glassy carbon with large graphene fragments. Shown in figure 6 is a snapshot of a small 
graphitic precursor molecule containing 3 fused-rings (two 6-membered rings and one 7-membered 
ring) formed after a 40 ps run at 3250 K. A commonly proposed reaction pathway [15] is that an open 
benzene ring and a closed benzene ring combine to form fused ring. However, our simulations illustrate 
that alkyl-like fragments from four different chains break, rearrange themselves and combine to form 



the fused ring fragment. Thus, it reveals another complex reaction pathway for graphitic precursor 
formation. This fused-ring structure is unstable at longer times, possibly due to the extremely high 
temperatures, and fragments into smaller carbon molecules. 

To identify a temperature zone where the graphitic precursors would be stable, detailed chain 
disintegration analysis is performed. In figure 7, the number of chains as a function of chain length 
(averaged over 4 different starting structures and a simulation time of 2 ps) is plotted on a semi-log scale 
for two different times - 20 ps (top panel) and 40 ps (bottom panel) and 3 different temperatures - 2750, 
3000 and 3250 K. At the start of the simulation run, the equilibrated structures have 16 chains with a 
carbon chain length of 56. At all temperatures after 20 ps, polymerization occurs between polymer 
chains resulting in longer chains of carbon chain length up to 120 at 3250 K and 190 at 3000 K. 
However, at longer times, all the long chains break to form smaller, shorter chains. This fragmentation 
into smaller chain occurs for other temperatures (2750 and 3000 K) at simulation times longer than 40 
ps. The number of C2 molecules (C2 stands for molecule which has two carbon atoms) is seen to 
increase. Included in the number of C2 molecules are a large amount of C 2 H radicals and C 2 gas 
molecules. This increase in C2 molecules with time is seen at all 3 temperatures. The fragmentation 
indicates that the temperatures are too high to polymerize and form large graphitic molecules. 

The results from figure 7 suggest that pyrolysis kinetics can be divided into two temperature zones. 
The high temperature (fragmentation) zone involves breaking of polymers into smaller molecules and 
the low temperature (polymerization) zone promotes bonding between polymers to form large 
connected graphitic network. In order to estimate a temperature range for these zones, in figure 8, the 
total number of C2 molecules formed after 40 ps is plotted as a function of temperature on a semi-log 
scale. The slope through the data points is extrapolated to intersect the x-axis (temperature-axis). The 
intersection occurs at 2050 K signifying a higher probability of graphitization at temperatures below 
2000 K. 

To capture large (stable) graphitic precursor formation a simulation was performed at 1750K. Since, 
no pyrolysis reactions are observed at 2000 K during the equilibration run for 250 ps, it can be expected 



that the time scale to observe any chemical reaction at 1750 K would go beyond the time scale 
accessible by MD simulations. Hence, to capture pyrolysis at 1750 K, a partially pyrolyzed structure 
from a 2750 K run after 40 ps was chosen as an input structure. In this structure more than 60% of the 
chains were intact. A reaction generally proceeds in 3 stages: initiation, propagation and termination. By 
using a partially pyrolyzed structure, the initiation stage is drastically reduced. In addition, after every 
1.25 ns, stable molecules such as H 2 0 and C0 2 are removed from the simulation cell. These stable gases 
interfere during the graphitization process preventing formation of large graphene precursors. At the end 
of a 12 ns simulation run two large graphitic structures are found, one containing 7 fused rings as shown 
in figure 9 and other containing 5 fused rings. The precursor starts with 2 fused rings and grows over the 
12 ns simulation run. This indicates that the reaction kinetics is different in the two temperature zones 
and it is necessary to perform the MD simulations at the temperature of experimental interest for 
accurate representation of the reaction behavior rather than merely increasing the MD simulation 
temperature to capture the reactions on the MD time-scale. 

3.2 Phenolic resin Composites 

Under this subsection, the effect of addition of CNT and carbon fibers on the initial stage of pyrolysis 
is discussed. 

A. CNT -Phenolic Resin Composite 

Pyrolysis simulations are performed at three different temperatures 2750, 3000 and 3250 K, on two 
composites containing CNT with (5,5) and (10,10) chirality. The CNT with two different chirality, is 
chosen to evaluate the effect of diameter on the pyrolysis simulation. In figure 10, the evolution of the 
primary product, water, is plotted as a function of time at three different temperatures, for CNT 
composites and compared with the pristine phenolic resin case. For all cases, the data is averaged over 4 
different starting structures to improve the statistical accuracy. The addition of fillers did not change the 
water production rate. The activation energy for water formation for the two types of CNT composites 
and the pristine case is same, 124 ± 20 kJ/mol. Thus, the chemical reactions in the initial stages of the 




pyrolysis are not affected by the presence of CNT, even at very high weight loading of 19% and 38% 
for CNT (5,5) and (10,10), respectively. 

Detailed analysis on the structure of CNT (5,5) after a simulation run of 40 ps for 4 different starting 
structures is performed. As shown in the various snapshots in figure 11, different starting structures 
yield a different number of chemically adsorbed molecules on the CNT surface. Hence, it is important 
for statistical purposes to perform 4 to 5 different simulation runs at each temperature of interest. Figure 
11(a), suggests that a CNT can play host to a wide range of molecules from a single hydrogen atoms to 
a 5-membered ring to large polymer chains. Even at 3000 K, the CNT is still stable after 40 ps and 
chemical adsorption of polymers. At 3250 K, in figure 11(c), CNT is distorted for 2 different starting 
structures, but it is found to be stable for other two starting structures after 40 ps. Separate simulations 
were performed on CNT without the phenolic resin at 3250 K to test the stability of CNT at this high 
temperature with ReaxFF potential. The CNT was stable at 3250 K during a 40 ps simulation run 
implying that the distortion of CNT is caused due to the physical and chemical interactions with the 
phenolic resin. 

B. Carbon fiber-Phenolic Resin Composite 

As mentioned in section 2, the carbon fiber is modeled as a block of graphite at the nanometer length 
scale. Pyrolysis simulation results are performed at 3250 K and averaged over 4 different starting 
structures. The number of water molecules formed is plotted against time and compared with the results 
from the pristine, CNT (5,5) and CNT (10,10) composite case in figure 12. A strong overlap of the 
water formation rate with no appreciable difference is seen in all cases. Similar overlap was seen at 
lower temperatures (2750 and 3000 K). Thus, it can be concluded that the presence of fillers such as 
CNT and carbon fibers do not play a catalytic role in the dehydration stage of the pyrolysis of the 
phenolic resin, even under high loading conditions. 

Even at 3250 K, for all the 4 structures studied, the graphite block is found to be stable during the 
simulation ran of 50 ps. At the end of the simulation run, small molecules such as methyl, ethyl group 
are bounded to the graphite layer. For the lowest pyrolysis simulation temperature studied at 2750 K no 
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molecules were bounded to graphite (in contrast to CNT, see figure 1 1(a)). For a simulation run at 3000 
K, similar to other two temperatures for the fiber-resin composite, no large molecules were found 
bounded to the graphite. Thus, in comparison to carbon fiber, a stronger covalent binding occurs in the 
case of CNT filler addition. 

The dehydration rate does not represent the entire kinetics of the phenolic resin pyrolysis and 
therefore it is unjustifiable to claim that the filler addition does not have any effect on the pyrolysis 
process. However, disintegration of polymer chains to form shorter molecules encompasses a large part 
of the overall pyrolysis process. Hence, we track the number of polymer chains with carbon chain length 
greater than 50 as a function of time for 4 different starting structures and at 3 different temperatures. 
The original carbon chain length of the phenolic resin in the starting structure is 56. The value of 50 
during the disintegration analysis is used for statistical purposes to prevent any chain end effects. At 
3000 and 3250 K, the disintegration rate for the pristine case (in blue) overlaps the CNT case (in green) 
and the graphite case (in red). At 2750 K, there is a slight discrepancy between the disintegration rates 
for simulation times greater than 20 ps. However, to confirm the effect of fillers, the simulations need to 
be performed in the low temperature zone (< 2000 K). Experimentally, it has been observed the addition 
of CNT produced cigar-like structures where the char strongly adheres to the CNT. 

4. Conclusion 

The primary objective of this paper was to model the effect of filler materials such as CNT and carbon 
fibers on the initial stage of the phenolic resin pyrolysis. Water was found to be the primary product for 
all the cases. There was no major effect of the filler addition on the water formation rate or the polymer 
disintegration rate. The evolution of graphitic precursors was studied rigorously. Two temperature zones 
were identified: the high temperature zone (>2000 K) involved fragmentation leading to non-graphitic 
structures and the low temperature zone (<2000 K) involved polymerization leading to predominantly 
graphitic structures. Formation of stable large (7 fused rings) graphitic precursors was captured at 1750 K 
during a simulation run of 12 ns. Thus, running the simulations at low temperatures was identified as the key 
to study the formation and growth of graphitic precursors. In the high temperature zone, no effect of fillers 



was observed on the graphitic precursor formation. To capture the cigar-like structures observed at 1000K 
experimentally for the CNT-resin composites, the simulations need to be performed in the low temperature 
reaction zone . 

The dehydration observed in our simulations is also found to be prominent in experiments during the 
curing stage of the resin which involves crosslinking of polymer chains to improve the mechanical 
properties of resin. However, the curing occurs at lower temperatures ~450 to 550 K. Jiang and 
coworkers showed that the simulation results for activation energy falls within the experimental range of 
74 to 199 kJ/mol for four mass-loss peaks [16]. The experiments were performed by Trick and Saliba on 
a cured sample of carbon-phenolic resin composite, where the pyrolysis kinetics were calculated from 
thermo-gravimetric analysis. It is difficult to compare our MD simulation results with such experiments 
because: (i) experiments are performed on a composite structure with carbon fibers as filler material; (ii) 
the phenolic resin is cured at 450 K in experiments, where most of the dehydration takes place resulting 
in crosslinking of polymer chains; (iii) for pyrolysis experiments, the heating rate is l°C/min in the 
temperature range of 573 to 1173 K, whereas simulations are performed with a step-change in 
temperature in the range of 2750 to 3250 K. Thus, in addition to a different starting structure, an 
additional assumption is needed that the reaction kinetics in the experimental temperature range of 573 
to 1173 K is same as the kinetics in the simulation temperature range of 2750 to 3250 K. According to 
our knowledge, there are no experiments to which we can directly compare our simulation results. 

The reactive MD simulations are definitely an important tool to provide insights on the involved 
reaction pathways and energy barriers. However, corresponding experiments are needed to validate 
these simulations . The experiments should include high heating rates to reach high temperatures (within 
a nanosecond time scale) such that the kinetics are monitored on the time scale accessible to MD 
simulations. State-of-the-art laser heating of polymer samples couple with spectroscopy is an option. In 
addition, the MD time scale needs to be increased from nano-second to micro-seconds in order to 
perform simulations at temperature of interest (1500-2000 K) without losing the atomistic information. 
This increase in MD time scale can be achieved by applying appropriate accelerated MD techniques 



[17]. Nevertheless, the results serve as a good model to provide insights on reaction mechanisms with 
caution while direct comparison with experiments. 
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Appendix A: Comparison of radial distribution function - ReaxFF vs. ab initio calculations 

Since ReaxFF is a relatively new force field, we tested its reliability for this specific application by 
comparing the radial distribution function (rdf) obtained from ReaxFF against results from ab initio 
based MD simulations. In particular, we performed simulations for a single phenolic chain with four 
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repeating units at T=300 K and confined in a periodic box of size 20 A on a side. The ab initio 
simulations were performed with the code VASP within the context of density functional theory (DFT). 
Note that unlike other implements of ab initio MD, VASP evolves the electronic and ionic degrees of 
freedom separately: forces obtained from DFT are converged first, and then used to advance the ionic 
positions. The GGA of Perdew and Wang was used with a cutoff of 300 eV and computations were 
done using the Projected Augmented Wave (PAW) potentials. Results are shown in Fig. 14. As can be 
seen the agreement is very good. Various bond lengths identified by peaks on the rdf agree very well 
although there is some variation in the peak heights. It should not be concluded that this represents a 
definitive test of ReaxFF. However, DFT is very good at describing bonded interactions, but is well 
known to have difficulties with weaker, long-ranged interactions that are common to polymer systems 
and are in principle included in ReaxFF. 

Appendix B: Water formation mechanisms 

Several dehydration mechanisms were observed during the different simulation runs. Similar to 
observations by Ref. 1, the dehydration mechanisms involved were: (i) between two -OH groups on 
phenols, (ii) a disassociated H and -OH group on a phenol, (iii) H from a benzene ring and -OH group 



from a phenol, and (iv) H from the methyl group (between two benzene rings) and -OH group on a 
phenol. For these four mechanisms, both intermolecular and intramolecular reactions took place. 

Shown in figure 15 is the temporal snapshot of a single polymer chain (for the pristine case) 
undergoing the fourth mechanism during intramolecular dehydration. The single polymer chain is 
isolated from the rest of the 15 polymer chains in the simulation cell for clarity. The simulation 
temperature is 2750 K and the simulation time is in the range of 9.025 to 9.2 ps. After 9.05 ps, the first 
transition state involves bonding between the -OH group from a phenol and hydrogen atom in the 
adjacent methyl group. At 9.075 ps, the -OH group breaks from the phenol, however, at 9.1 ps, the 
water molecule reattaches to the benzene ring. After 9.2 ps, the water molecule breaks free of the 
polymer and diffuses away from the parent phenolic resin chain. During this period, the polymer chain 
is intact without undergoing any substantial disintegration (except thermal dissociation of couple of 
hydrogen atoms). Thus, fig. 15 provides an atomistic picture of the transition states involved during a 
particular water formation mechanism for a thermally stable polymer chain. This process can be adopted 
to develop reaction pathways and extract associated chemical kinetic information to build highly 
detailed atomistically-informed chemical kinetic models. 
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FIGURE CAPTIONS 


Figure 1. Top Panel - Snapshot of a single phenol formaldehyde polymer chain with 8 repeating units 
and terminated with a methyl group on one end; color code: Carbon - green, Hydrogen - blue, Oxygen 
- brown. Bottom Panel - Snapshot of the simulation cell containing 16 phenol formaldehyde polymer 
chains equilibrated at 300 K after annealing. 

Figure 2. (a) Snapshot of equilibrated resin composite structure with CNT (5,5) at the center of the 
simulation cell; (b) Snapshot of equilibrated resin composite structure with a graphite block (containing 
3 graphene layers in A-B-A stacking sequence) at the center of the simulation cell. The color code for 
both snapshot is: Carbon - green, Hydrogen - blue, Oxygen - brown. 

Figure 3. Total radial distribution function (rdf) at 300K. The results match qualitatively with pcff 
forcefield used in Ref. 14. 

Figure 4. (a) Number of water molecules (averaged over 4 different simulation runs) plotted against 
time for 5 different temperatures - 2750 K (red), 2875 K (blue), 3000 K (green), 3125 K (black) and 
3250 K (magenta), (b) Arrhenius plot of the logarithm of the water formation rate (k) versus 
temperature gives an activation energy of 124 ± 20kJ/mol. This value matches well with Ref. 1. 

Figure 5. Number of water molecules produced as a function of time for 2750 K (blue), 3000 K (green) 
and 3250 K (red). For all temperatures at long simulation times, the number of water molecules 
produced asymptote to a value near 100. 

Figure 6. Snapshot of a small graphitic precursor molecule containing 3 fused rings out of which one is 
a 7-membered ring. The molecule was found after a 40 ps run at 3250 K. The atoms that form this 
molecule come from 4 different chains. Color code: Carbon - green, Hydrogen - blue, Oxygen - brown. 



Figure 7. Average number of chains as a function of carbon chain length on a semi-log scale at 3 
different temperatures, 2750, 3000 and 3250 K and at two different times, 20 ps (top panel) and 40 ps 
(bottom panel). Disintegration of polymer chains to form small molecules is prominent at all 
temperatures. 

Figure 8. The number of C2 molecules formed after 40 ps (averaged over 4 different starting structures) 
is plotted as a function of temperature on a semi-log scale. On extrapolating the slope, it intersects on 
the x-axis at 2050 K, suggesting the existence of two temperature zones, disintegration and 
polymerization. 

Figure 9. Snapshot of a large graphitic precursor structure formed after a simulation run time of 12 ns at 
1750 K. Lower temperatures (< 2000 K) are necessary to capture the formation and growth of graphitic 
precursors. Color code: Carbon - green, Hydrogen - blue, Oxygen - brown. 

Figure 10. Number of water molecules produced (averaged over 4 different starting structures) as a 
function of time for 3 different cases: pristine, CNT (5,5) composite and CNT (10,10) composite at 
2750, 3000 and 3250 K. The pristine case is represented in blue, CNT (5,5) composite in green and 
CNT (10,10) in red. No effect of CNT filler on water production rate is seen. 

Figure 11. Snapshot of CNT (5,5) after a simulation run of 40 ps for 4 different starting structures at 3 
different temperatures, 2750, 3000 and 3250 K. At 3250 K, the CNT is found to be distorted for 2 
different starting structures. The distortion results from physical and chemical interactions with the 
resin. Color code: Carbon - green, Hydrogen - blue, Oxygen - brown. 

Figure 12. Comparison of the number of water molecules formed (averaged over 4 different starting 
structures) against time at 3250 K for the pristine, CNt (5,5), CNT (10,10) and graphite case. No 
substantial effect of filler addition is seen on water formation rate. 
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Figure 13. Number of phenolic resin chains with carbon chain length greater than 50 (averaged over 4 
different starting structures) is plotted as a function of time for 3 different cases: pristine, CNT (5,5) 
composite and graphite composite at 2750, 3000 and 3250 K. The pristine case is represented in blue, 
CNT (5,5) composite in green and graphite composite in red. This shows the disintegration rate of the 
polymer chains. 

Figure 14. RDF for a short, four unit phenolic chain comparing the force field, ReaxFF versus ab initio 
molecular dynamics at 300 K. 

Figure 15. Temporal snapshot of intramolecular dehydration mechanism between -OH group on phenol 
and hydrogen from adjacent methyl group. The polymer chain stays intact during the reaction period. 
Color code: Carbon - gray, Hydrogen - white, Oxygen - red (for atoms involved in reactions, Oxygen - 
blue, Hydrogen - magenta). The region of interest is highlighted by a square white box. 
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